#include <numeric>
#include <kuai/misc/error.hpp>
#include <kuai/sim/reporter.hpp>

namespace kuai {

	SpaceReporter::SpaceReporter() 
	{
		pout = NULL;
		dumpLevel = DUMP_DETAIL;
	};


	void SpaceReporter::set(const String& key, const String& data) {
		if (boost::algorithm::iequals(key, "logfile")) {
			filename = FileName(data);
			tryOpen(fout, filename);
			pout = &fout;
		}
		else if (boost::algorithm::iequals(key, "level")) {
			if (boost::algorithm::iequals(data, "off") || boost::algorithm::iequals(data, "nothing")) {
				dumpLevel=DUMP_NOTHING;
			}
			else if (boost::algorithm::iequals(data, "summer")) {
				dumpLevel=DUMP_SUMMER;
			}
			else if (boost::algorithm::iequals(data, "detail")) {
				dumpLevel=DUMP_DETAIL;
			}
			else if (boost::algorithm::iequals(data, "all") || boost::algorithm::iequals(data, "everything")) {
				dumpLevel=DUMP_EVERYTHING;
			}
			else {
				throw error("Unknown dump level %1%", data);
			}
		}
		else {
			throw UNDEFINED_MEMBER(key);
		}
	};
	
	String SpaceReporter::get(const String& key) const {
		if (boost::algorithm::iequals(key, "logfile")) {
			return str(filename);
		}
		else if (boost::algorithm::iequals(key, "level")) {
			switch (dumpLevel) {
			case DUMP_NOTHING:
				return "nothing";
			case DUMP_SUMMER:
				return "summer";
			case DUMP_DETAIL:
				return "detail";
			case DUMP_EVERYTHING:
				return "everything";
			default:
				assert (false);
				break;
			}
		}
		
		throw UNDEFINED_MEMBER(key);
	}

	IMPLEMENT_WORKER_CREATOR(CountAngleDistribution, angle_distribution)


	CountAngleDistribution::CountAngleDistribution() {
		axis = XYZ(0, 0, 1);
		pbc = NULL;
		slice = deg2rad(1);
	}

	void CountAngleDistribution::set(const String& key, const String& data) {
		if (key == "atoms") {
			std::istringstream stream(data);
			Index a1, a2;
			while (stream >> a1 >> a2) {
				pairs.push_back(std::make_pair(a1-1, a2-1));
			}
		}
		else if (key == "slice") {
			slice = deg2rad(lexical_cast<RealNumber>(data));
		}
		else if (key == "axis") {
			axis = lexical_cast<XYZ>(data);
		}
		else {
			SpaceReporter::set(key, data);
		}
	};
	
	String CountAngleDistribution::get(const String& key) const {
		if (key == "atoms") {
			std::ostringstream stream;
			for (size_t i = 0; i < pairs.size(); ++i) {
				stream << pairs[i].first << " " << pairs[i].second << std::endl;
			}
			return stream.str();
		}
		else if (key == "slice") {
			return str(rad2deg(slice));
		}
		else if (key == "axis") {
			return str(axis);
		}
		else {
			return SpaceReporter::get(key);	
		}
	}

	void CountAngleDistribution::start(std::ostream& stdlog, std::ostream& errlog, PhaseSpace& space, EnergySet& efunc, Index tid, Index nThreads) {
		if (tid == 0) {
			if (pout == NULL) {
				pout = &stdlog;
			}

			PhaseSpace::const_iterator it = space.find(MOL_INFO_ATOM_COORD_3D);
			coords = static_cast<const PropertyXYZ*>(it->second.get())->c_ptr();
			if (it == space.end()) {
				throw UNDEFINED_MEMBER("ATOM_COORD");
			}

			it = space.find(MOL_INFO_PBC);
			if (it != space.end()) {
				pbc = static_cast<const PropertyPBC*>(it->second.get())->c_ptr();
			}

			if (pairs.empty()) {
				throw UNDEFINED_MEMBER("atoms");
			}
			else {
				size_t nAtoms = space.countAtoms();
				for (size_t i = 0; i < pairs.size(); ++i) {
					if (pairs[i].first >= nAtoms || pairs[i].second >= nAtoms) {
						throw error("Atom Pair %1%, %2% is out of range!", pairs[i].first, pairs[i].second);
					}
				}
			}

			count.resize(size_t(3.142 / slice), 0);

			if (dumpLevel >= DUMP_DETAIL) {
				*pout << 
					"============================================================\n" 
					"       frame |  atom1  atom2 | degree | column\n"
					"------------------------------------------------------------" 
					<< std::endl;
			}
		}
	}
	
	void CountAngleDistribution::act(PhaseSpace& space, size_t iStep, EnergySet& efunc, Index tid, Index nThreads) {
		if (tid == 0) {
			for (size_t i = 0; i < pairs.size(); ++i) {
				XYZ delta = coords[pairs[i].first] - coords[pairs[i].second];
				if (pbc) {
					pbc->update(delta);
				}
				RealNumber theta = angle(delta, axis);
				int index = int(theta/slice);
				if (index < 0) {
					index = 0;
				}
				else if (index >= int(count.size())) {
					index = count.size() - 1;
				}
				++count[index];
				if (dumpLevel >= DUMP_DETAIL) {
					char buf[128];
					sprintf(buf, "%12u | %6d %6d | %6.2f | %6d", 
						iStep, pairs[i].first+1, pairs[i].second+1, rad2deg(theta), index);
					*pout << buf << std::endl;
				}
			}
		}	
	}

	void CountAngleDistribution::finish(std::ostream& stdlog, std::ostream& errlog, PhaseSpace& space, EnergySet& efunc, Index tid, Index nThreads) {
		if (tid == 0)
		{
			if (dumpLevel >= DUMP_DETAIL) {
				*pout << 
					"============================================================\n\n\n\n" 
					<< std::endl;
			}
			if (dumpLevel >= DUMP_SUMMER) {
				char buf[128];
				*pout << 
					"                        Angle Distribution\n"
					"============================================================\n" 
					"       | a1(deg)  a2(deg)|        count |    %\n"
					"------------------------------------------------------------\n";

				size_t nTotal = std::accumulate(count.begin(), count.end(), 0);

				for (size_t i = 0; i < count.size(); ++i) {
					RealNumber a1 = rad2deg(slice*i);
					RealNumber a2 = rad2deg(slice*(i+1));
					RealNumber percent = count[i] * 100.0f / nTotal;
					sprintf(buf, "%6u | %6.2f ~ %6.2f | %12u | %6.2f%%\n", i+1, a1, a2, count[i], percent);
					*pout << buf;
				}

				*pout << 
					"============================================================\n\n\n\n" 
					<< std::endl;
			}
		}
	}
	

	WORKER_HELP_FUNCTION(CountAngleDistribution, output) {
		// TODO:
	}


}